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1 . INTRODUCTION 


1 .1 Purpose of this Report 

Two design objectives for Stirling engines are to keep mechanical 
losses low and to achieve good heat transfer from and to the working fluid. 
These objectives conflict because heat exchangers with higher heat transfer 
effectiveness generally have higher flow losses. It is therefore valuable 
for the designer to know the trade-off between pressure drop and heat 
transfer. This Is particularly important when high performance is required 
as :'r the Automotive Stirling Engine and the Space Power Demonstrator 
Engine. 

Computer modelling is used to evaluate options in Stirling engine 
design. Measured Stirling engine efficiencies and power output are often 
lower than predictions would indicate, however. Transient pressure 
variations in the Stirling engine working spaces are also often not 
predicted correctly (Rix 1 98*4 ) . Inadequate modelling of fluid mechanics and 
heat transfer in the heat exchangers (usually by using steady-state 
correlations for pressure drop and heat transfer) may partially account for 
this discrepancy. Improving this model will require more physically correct 
expressions for the effect of oscillation. It is hoped that this report 
will make a step in that direction by improving the understanding of fluid 
mechanics and heat transfer in Stirling engines. It summarizes, extends and 
discusses research found in the literature and draws tentative conclusions. 

1 .2 Outline 

The oscillating flow effect is isolated as the focus of this report. 
Similarity parameters to characterize fluid mechanics and heat transfer in 



Stirling engine heat exchangers are proposed. A literature review and some 
analyses address fluid mechanics and heat transfer in oscillating flow in 
pipes under conditions which are similar to conditions in heaters and 
coolers. Literature which deals with flow in porous media, similar to 
regenerators, is also reviewed. The operating characteristics of 11 
Stirling engines are described in terms of similarity parameters and 
tentative conclusions are drawn about the fluid mechanics and heat transfer 
conditions in Stirling engine heat exchangers. Open questions to be 
answered by future research are identified. Some of the material covered in 
this report is summarized in Seume and Simon (1986). 
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2. SIMILARITY PARAMETERS 


Similarity parameters are required to concisely and generally describe 
fluid mechanics and heat transfer in physically similar systems. Their 
proper choice is crucial for the development of an experimental program and 
the interpretation of experimental data. 

Section 2.1 discusses some characteristics of the flow in Stirling 
engine heat exchangers that are not described by similarity parameters — in 
particular, simplifying assumptions to isolate the oscillating flow effect. 
Section 2.2 presents physical arguments for the chosen set of similarity 
parameters and Section 2.3 derives the similarity parameters by normalizing 
the governing equations. 

2.1 Isolating the Oscillating Flow Effect 

The effect of oscillation on fluid mechanics and heat transfer has been 
the subject of several studies related to Stirling engines (e.g. , Kim 1970, 
Organ 1975, Miyabe et al. 1982, Chen and Griffin 1983, Hwang and Dybbs 1980 
and 1983, Taylor and Aghili 1984, Dijkstra 1984, Jones 1985, Rice et al. 
1985). In review of these studies, it became clear that the effect of 
oscillation on pressure drop, viscous dissipation and heat transfer must 
first be studied isolated from the remainder of the processes occurring in 
the engine, in particular isolated from: 

- compression and expansion of the working fluid 

- non-harmonic fluid motion 

- high temperature gradients. 

Some discussion about each of these effects in the Stirling engine follows. 
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Compression and expansion of the working fluid . The working fluid in 
Stirling engines moves from the compression space toward the expansion space 
and back during each cycle. Simultaneously, it undergoes compression and 
expansion (Figure 2-1). Therefore, the mass of fluid in the heater in 
Figure 2-1 varies with time and consequently the mass flow rates at the ends 
of the heater differ throughout the cycle (full and dashed line). 

Appendix A provides a survey of the cycle variations of pressure and 
Reynolds numbers in Stirling engine heat exchangers using the isothermal and 
Schmidt analyses (c.f. Section 6.1). The engines and their operating points 
are described in Section 6.2. Figure 2-1 shows that the pressure change is 
fastest when the fluid velocity in the heater is near zero. In the cooler 
and in the regenerator, however, velocities are high when the pressure is 
changing most rapidly. 

During compression, the gas temperature rises. This decreases the 
temperature difference between the bulk of the gas and the heater wall and 
thereby may reduce the heat transfer from the heater pipe to the gas. In 
the cooler, compression increases the temperature difference between the 
bulk fluid temperature and the wall temperature. The heat transfer may, 
however, be out of phase with the bulk-to-wall temperature difference. This 
was shown by Faulkner and Smith (1984) in a study related to heat transfer 
in the cylinders of reciprocating machinery. In the heater, fluid 
velocities are low during expansion and compression. Consequently, 
convective heat transfer is not expected to be as strong as in the cooler 
and regenerator during compression and expansion. 

In the cooler and regenerator, temperature in the gas increases 
simultaneously with convective heat transfer during the cold blow, i.e., 
when gas moves from the compression space towards the expansion space. 
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Figure 2-1: Variation of pressure and Reynolds nunber 

in the heater of a GPU 3, isothermal analysis. 

Flow toward the cold end is positive; 

full line *= hot end, dashed line = cold end. 


5 




During the hot blow, the gas expands as it flows through the cooler and 
regenerator . 

Dijkstra (198H, p. 1886) proposed to model the compression of the 
working gas as bulk heating of a fluid flowing through heat exchangers. He 
pointed out that the convective heat transfer coefficient in the case of 
bulk heating is very different from that for specified wall temperature or 
specified heat flux. The differences are particularly big in the entrance 
region of tubes and in the turbulent flow regime. It is, however, not clear 
whether bulk heating can be used to model temperature changes in the working 
gas due to compression and expansion. 

This report focuses on the effect of the oscillation of the working 
fluid on fluid mechanics and heat transfer in the heat exchangers. To 
isolate this effect, compression and expansion of the working fluid are 
neglected in most of the discussion below. 

Mean velocity variation . In Figure 2-2, the mean velocity is 
normalized with the speed of sound to form a Mach number; according to the 
assumptions implicit in the isothermal analysis, the speed of sound is 
uniform and constant in each heat exchanger. Figure 2-2 shows that the mean 
velocity variation in a typical Stirling engine is roughly sinusoidal and 
very similar to the Reynolds number variation shown in Figure 2-1 . 

Velocity and Reynolds number variation deviate from a sine function in 

that: 

(1) the hot blow period, during which the flow is positive (i.e., 
towards the compression space), is shorter than the cold blow 
period. 

(2) The maximum velocity during the hot blow period is greater than 
during the cold blow period. 
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The differences of the mass flow variation (Figure 2-1) and the mean 
velocity variation (Figure 2-2) from harmonic motion will be neglected in 
the discussion below. 

Temperature dependent properties . The molecular viscosity and thermal 
conductivity of gases increase with temperature. Due to the temperature 
gradients in the heat exchangers, the gas properties vary in the radial and 
axj.al directions. 

In particular, the viscosity near the heater wall is expected to be 
higher than in the core; the viscosity of the ga3 is lower near the wall of 
the cooler than in the core of the flow. Therefore calculations based on 
the assumption of constant properties in this report should only be 
considered approximate. Some authors consider density variations (buoyancy 
effects) in turbulent forced convection with large density gradients 
important. For a review of this subject, see Petukhov et al. (1982). 

2.2 Physical Arguments for the Choice of Similarity Parameters 

This section describes qualitatively some of the transient phenomena 
expected in Stirling engine heat exchangers. Similarity parameters 
characterizing each of these phenomena are proposed. 

Velocity and acceleration . In steady flow, the friction factor (a 
non-dimensional viscous pressure drop) is a function of Reynolds number (a 
non-dimensional mass flux). Oscillating flow experiments, e.g. those by 
Taylor and Aghili ( 1 984 ) showed that the functional relationship between 
friction factor and Reynolds number is different for oscillating flow than 
for steady unidirectional flow. 

In the limit of very slow oscillations, the oscillating flow 
relationships should approach that of steady flow. This suggests that the 
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friction factor in oscillating flow is also a function of a dimensionless 
frequency or dimensionless acceleration. Following Chen and Griffin (1983) > 
this dimensionless frequency is chosen as: 


Re 

u 


a) d 
H v 


2 


In this report, the Reynolds number is based on the maximum velocity in the 
cycle : 


R e max 


u m»max d 

v 


Fluid displacement . Walker (1962, 1963 and 1980, p. 130) shows that 
the fluid displacement in Stirling engines may be so small that some working 
fluid merely moves back and forth within the heat exchanger without exiting. 
If the gas moves as a plug and if there is no axial mixing, this gas 
contributes to the axial heating only by absorbing thermal energy from the 
wall or adjacent fluid at the high temperature end of the heat exchanger and 
releasing it at the low temperature end. Only the gas that enters and 
leaves the heat exchanger transfers heat directly to the volume in the 
cylinders. Organ (1975) argued that the gas that stays within the heat 
exchanger may be considered an undesirable addition to the dead space in the 
engine. The amplitude ratio, Ap, is used to describe fluid displacement in 
heat exchangers. It is the fluid displacement during half a cycle divided 
by the tube length, computed by assuming that the fluid moves as a plug at a 
mean velocity, u m . This is the reciprocal of A introduced by Organ (1975). 
If Ap « 1, most fluid oscillates within the tube without exiting; if Ar » 1, 
the fluid traverses quickly through the tube, residing in the upstream and 
downstream spaces during most of the cycle. 
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Entrance and exit effects . The pressure drop at the entrance of a heat 
exchanger has two components. The reversible pressure drop is the result of 
an Increase in kinetic energy by expending flow work (pressure drop) in the 
contraction upon entrance to a heat exchanger. In the case of a 
compressible working gas, density will decrease as the gas is expanded 
(approximately adiabatically) in this contraction. The irreversible part of 
the entrance pressure drop is due to viscous dissipation in the shear flow 
of the vena contracta. It is typically small in steady flows. At the exit, 
pressure may be recovered as the gas velocity decreases. The pressure rises 
and kinetic energy is lost irreversibly by viscous dissipation in the shear 
layers of the separated flow upon exit from the heat exchanger into the 
downstream chamber. 

Since the magnitude of entrance and exit losses depends on the area 
ratios of the contraction and expansion of the duct, these area ratios (or 
the corresponding diameter ratios, (d/D)) are added to the Reynolds number 
as similarity parameters for capturing entrance and exit losses. 

Pipe curvature . The pipes in Stirling engine heaters are often curved. 
It is known from steady unidirectional flow that secondary flows develop in 
curved pipes as centrifugal forces act on the fluid. The Dean number, 

d ® " Re /sr* 

where R c is the radius of curvature of the pipe centerline, is commonly used 
as a similarity parameter for this effect (Berger et al. 1983 ). 

Developing flow . In steady flow, the velocity profile in a pipe 
changes in the flow direction until fully developed conditions are reached. 
The entrance length scales on the diameter of the pipe (or the hydraulic 
diameter of the duct). The friction factor in the entrance (developing) 
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region is higher than in the fully developed region. Therefore, the average 
friction factor for developing duct flow is a function of the 
length/diameter ratio, i/d, as well as the Reynolds number. 

Flow in regenerators . Flow in the complex geometries of Stirling 
engine regenerators is commonly described in terms of a Reynolds number 
based on the hydraulic diameter and the matrix porosity. Other choices of 
similarity parameters are discussed in Section 5.1. 

Compressibility of the working gas . The compressibility of the working 
gas affects the entrance and exit losses of the heat exchangers as described 
above and the propagation of pressure changes throughout the engine. If the 
gas velocity in a heat exchanger is of the same order as the speed at which 
pressure changes propagate (approx, the speed of sound), pressure (and, 
therefore, density) variations in the engine due to compressibility will be 
significant. As the ratio of gas velocity to speed of sound (i.e. the Mach 
number) approaches unity, shock waves form and the flow becomes choked. The 
Mach number, based upon the highest velocity in the engine cycle, is 
therefore chosen as a measure of compressibility effects. 

Temperature transients . The walls of the heater and cooler and the 
matrix of the regenerator undergo temperature transients as they are 
subjected to different gas temperatures during hot and cold blows. The 
transient response, particularly of the regenerator matrix, affects the 
engine performance (Walker 1980, p. 140). This response is described by 
dimensionless parameters from regenerator theory, e.g. dimensionless 
regenerator length and dimensionless blow period, as discussed in Section 
5.3. Also, the transient response of the heater and cooler walls may 
contribute to enhanced axial heat transfer as discussed in Section il.l. 
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2.3 Derivation of Similarity Parameters 


Momentum equation . Similarity parameters, introduced in Section 2.1, 
can be derived by normalizing the momentum and the energy equations. 
Neglecting gravity, the momentum (Navier-Stokes) equation is: 

. 3.*; . - *E . vV 2j 

at p 

A normalized form suitable for pipe flow is derived in Appendix B: 

Re Jh! + Re max +*. v «+* - - Re max ^R* + v V*u* 

« at* 2 2 p 

Here the superscript * denotes a normalized quantity. The dimensionless 
frequency, Re u , is the coefficient of the temporal acceleration (or 
unsteady) term, whereas the maximum Reynolds number , R^max* is the 
coefficient of the spatial acceleration and the pressure gradient term. 

The viscous term has no coefficient. This does not imply, however, 

that Re max is the ratio of the u-Vu (steady inertia) term to the vV u 

(viscous) term. That would only be true if both non-dimensional terms were 
of the same order, i.e., if 

(Xu'-Vi?) - 0(v*V*’u*) 

This is not necessarily the case. In a fully developed laminar pipe flow, 
for example, the u*Vu term vanishes but the Reynolds number is still 

non-zero. 

The dimensionless frequency, Re^, has also been called the kinetic 
Reynolds number (White 1974, p. 144) or the Valensi number (Park and Baird 
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1970). It Is a multiple of the Stokes number (Grassmann and Tuma 1979) and 
it is the square of the Womersley parameter. 

A slightly rearranged version of the normalized momentum equation is: 


Str 3u* _-*•* V*p* 2 

— -r-r + u *Vu Jr + =- v V 2 u 

2 3t p Re m ax 

where Str - ud/u m|inax is the Strouhal number. Note that the definition of 
the Strouhal number may vary in the literature (Telionis 1981, p. 88), and 
that Str is a multiple of the Strouhal number used by Dijkstra (1984). 

Geometric similarity . The choice of length scales for the 
normalization of the momentum equation is arbitrary. To maintain 
similarity, however, all dimensions must scale on the same length; in this 
case the pipe radius or the diameter d was chosen. Therefore, (fc/d) is the 
appropriate similarity variable describing the pipe length. The relative 
amplitude of fluid displacement is not an independent similarity parameter. 
For sinusoidal fluid motion, it is shown in Appendix B to be: 


lB-53 


2 d Re 


4ll 


Description of geometries other than straight pipes, such as the 
various regenerator types, curved pipes, or ducts of rectangular 
cross-section, require additional descriptors to complement the similarity 
parameters described in this section, e.g. the Dean number. 

Energy equation . Heat transfer is governed by the energy equation: 


pc p (f ♦ S-vn 


i£ 

at 


+ u*Vp + V*(kVT) 


+ 


* 
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where 4> is the dissipation function (Appendix B). A normalized form for 
pipe flow is: 


p*c p * (2 Re^ ■|^» + Re max u*-7*T*) 


2 Re^ Ec ||* ♦ Ec Re max u*-V*p 


+ vMkW) + Ec 

p r 


where Pr - C PQ 
K o 


Ec 


Urn. 


mimax 


c po (Th”T c ) 


Other choices of similarity variables . Other forms of the normalized 
governing equations are obtained if the variables are normalized differently. 
If the temperature is normalized as T* - T/T 0 and we assume that the fluid 
is an ideal gas with constant specific heats, the Eckert number, Ec, can be 
replaced by (Y - 1) M max 2 . 

The Mach number can also enter the normalized momentum equation through 
the equation of state if pressure is normalized with a reference pressure: 


_ P_ 
Po 


Po “ PoRTo 


Pq ( u m,max)2 
V Mmax 


Then the Mach number appears in the pressure term: 


Rp 


3u* 


_ , „# * 

Re maY 1 V p 

~o ym r + v V 2 u 

2 >Mmax P 


These examples show that the choice of similarity parameters is 
somewhat arbitrary. The appearance of a dimensionless parameter in a 
normalized equation alone does not imply that it is a relevant similarity 
parameter. Meaningful similarity parameters are obtained when all 
normalized terms in an equation are of order 1 , e.g. in the momentum 
equation: 


14 




expected. 
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3. FLUID MECHANICS OF OSCILLATING PIPE FLOW 


In unsteady flows, pressure and shear forces in the fluid are not in 
equilibrium; instead there is a local balance of inertia, pressure and shear 
forces. Therefore, the fluid mechanics of oscillating flows are different 

from those of steady flow. 

Figure 3-1 presents an integral momentum balance on a cylindrical 
control volume in fully developed flow (u-Vu - 0). The oscillating flow 

case shows the additional term associated with temporal acceleration (9u/9t, 
where u is the cross-sectional mean velocity). 
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(a) pressure in balance with wall shear stress 


w 


VD 2 p( 


x ) f A D‘ : 'ax8 u/3 1 /*" 

_i rJ ~ d 


&ttD p(x+Ax) 


oscillating flow 


(b) pressure in balance with wall shear stress + fluid inertia 


Figure 3 - 1 : Force balance in fully developed pipe flow. 
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3.1 Analysis of Laminar Flow In a Straight Pipe 

Measurements by Richardson and Tyler (1929) first indicated an 
oscillating flow effect. They found a maximum velocity near the wall for 
oscillating pipe flow. Sexl (1930), Womersley (1955) and Uchida (1956) have 
since shown this by analysis of both sinusoidal and non-sinusoidal 
oscillating flows. Similar analyses were performed by Drake (1965) for 
rectangular channels and by Gedeon (1986) for flow between parallel plates. 
Zielke (1968) and Trikha (1975) used Laplace transforms to calculate 
pressure drop in unsteady laminar flow. 

Figure 3-2, based upon the Uchida (1956) analysis, shows velocity 
distributions for fully-developed flow at several times during the cycle and 
for several values of Re w . All profiles assume the same magnitude of 
pressure gradient oscillation, and the local velocities are normalized with 
the maximum mean velocity that would occur with steady laminar flow 
(responding to the cycle-maximum pressure gradient). Velocity distributions 
similar to these were measured by Shizgal et al. (1965, p. 99) and Edwards 
and Wilkenson (1971, p. 87). An analysis similar to Uchida's but for flow 
between parallel flat plates was made by Gedeon (1986). Figure 3~2 shows 
parabolic distributions for slowly oscillating pipe flow (Re^ - 1). With 
higher Re u , however, the velocity amplitude decreases and, during parts of 
the cycle, the near-wall flow direction is opposite to the mean flow 
direction (Re,,, - 30). A velocity maximum develops near the wall. At the 
maximum frequency (Re u - 1000), this free shear layer becomes narrower and 
moves closer to the wall. A uniform velocity core exists and velocity 
gradients become concentrated in a Stokes layer near the wall. 
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Figure 3-2: Velocity profiles in oscillating, laminar, fully developed pipe flow 


The presence of a large velocity gradient between the maximum velocity and 
the core (the free shear layer) may play an important role in 
laminar-to-turbulent transition. Further velocity profiles are shown in 
Appendix C. 

To discuss the effect of flow oscillation on shear stress and pressure 
drop, augmentation factors, a T and ap, and lead angles, X T and Ap, are 
introduced as: 


u m “ 

u m, 

max cos 

( U)t ) 


*W “ 

“t 

8|i/d u m| 

max cos 

(wt + A t ) 

Ap - 

°P 

32 pl/d* 

u m,raax 

cos (wt + Xp) 


Note that for steady flow, a T and a p are unity and X T and X p are zero. 

These quantities, taken from the Uchida (1956) analysis, are plotted versus 
Re w on Figures 3-3 and 3-4. Note (Fig. 3-3) that for the Re^ range 
chosen, the wall shear stress is enhanced by a factor of 8 over that of 
unidirectional flow, and that the pressure change is enhanced by a factor of 
130 . 

The phase relationships given in Figure 3“ 4 indicate that, in the limit 
of high Re u ,the shear stress leads the mean velocity by 45° and the 
pressure change leads the mean velocity by 90°. The near-wall fluid 
responds more quickly to the imposed pressure gradient than does the higher 
inertia core flow. Clearly, neither the shear stress nor the pressure drop 
should be computed from the pressure gradient by assuming quasi-steady flow. 

In oscillating flow, loss of engine work due to viscous dissipation 
cannot be calculated by multiplying mean velocity and wall shear stress as 
is done in steady flow. Instead, the viscous power dissipation in a pipe 
section of length i, must be computed as: 
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Figure 3-3: Coefficients of amplitude of pressure drop 

and wall shear stress 
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t o o r d/2 / 2 m 

\ m 2 * 1 K y( aF } r dr 

The relative increases (due to flow oscillation) of shear (irreversible) 
work and total pumping work (reversible and irreversible) over those for 
steady flow can be expressed by the coefficients of excess work, and e p , 
respectively. 

i 

[J 2 J rut) d(tat) ] 

[/ o* I ,(“t) d(»t)] Re ^ 

0) 

£il ~ ( &f >u . dUt) l 

CP ‘ t/ V - «*>]*<, 

< b) 

These quantities are plotted versus Re u on Figure 3-5. Note that over the 
Re w range shown, the' irreversible work increases by a factor of six, 
whereas the pumping work increases by a factor of 125. The difference 
between the pumping and shear work values represents reversible work. In 
the Stirling engine most of this is eventually lost due to dissipation upon 
entry to or exit from the heat exchanger channels. Some of this reversible 
* work may be recovered as gas exiting the channels impinges upon and does 

work on a piston surface. 
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Figure 3-5: Relative increase of pumping work due to flow oscillation 
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3.2 Laminar Oscillating Flow In Curved Pipes 


The pipes in many Stirling engine heaters are bent. If the radius of 
curvature of a pipe is small, secondary flows are strong. This is true for 
steady and for oscillating flows. Due to the pipe curvature, centrifugal 

forces act on the fluid; therefore the spatial acceleration term u*Vu is 

non-zero while it is zero in fully developed laminar flow in a straight 
pipe. 

In the normalized momentum equation, the Dean number replaces the 
Reynolds number. The unsteadiness is again characterized by Re w . 

Telionis (1981 , p. 1 83 ) gives a different but equivalent set of similarity 
parameters. Figure 3-6 shows a schematic of oscillating flow in a curved 
pipe. In each half of the pipe there are two counter-rotating secondary 
flows: one in the potential core and one in the Stokes layer. 
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a - d/2 


r o “ R c 


Figure 3-6: Oscillating flow in a curved pipe 

(from Telionis 1981, p. 182) 


Yamane et al. (1985) calculated these streamlines and axial velocity 
profiles. The latter were confirmed experimentally by Sudou et al. (1985). 
Streamlines and velocity profiles for laminar, oscillating flow in 
rectangular curved ducts were calculated by Sumida and Sudou (1985). 

Yamane et al. (1985) calculated and measured the increase of wall shear 
stress and pressure drop in a curved pipe over that in a straight pipe. 

They presented these ratios as a function of Dean number and Womersley 
parameter. 
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3.3 Non-sinusoidal and Free Oscillations 


Section 3.1 discussed sinusoidal oscillations in a straight pipe. 

Uchida (1956) used a Fourier series representation to extend this analysis 
to general, non-sinusoidal periodic laminar flows. This analysis may be 
useful since oscillations in Stirling engines are generally not sinusoidal. 
Chan and Baird (1974) studied forced oscillations of liquid columns in a 
range that is of interest for the design of liquid piston engines (West 
1983, pp. 54-59, 131). Free (damped) oscillations of liquid columns in 
U-tubes were studied by Valensi (1947), Park and Baird (1970), and Iguchi et 
al. (1982). It is not clear to what extent their results regarding fluid 
friction and the transition from laminar to turbulent flow (see Section 3.4) 
can be applied to the forced oscillations of liquid columns occurring in 
liquid piston engines. 

The effect of pipe curvature on the flow was neglected in the studies 
referred to above and in the transition studies discussed in the next 
section. The authors assumed that the curvature effect is negligible 
because the pipe radius of curvature is small compared to the pipe radius. 
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3.4 Transition from Laminar to Turbulent Flow 

The transition process . Transition in unidirectional steady flow is 
known to be sensitive to bulk mean velocity and acceleration. Lower 
velocities and acceleration stabilize, whereas higher velocities and 
deceleration destabilize. Since in oscillating flows both velocity and 
acceleration vary, it is expected that flow patterns may change from 
laminar-like to turbulent-like throughout the cycle under near-critical 
conditions. Hino et al. (1976) probed an oscillating pipe flow with a 
single hot-wire and presented the traces of absolute value of the velocity 
shown in Figure 3-7 . The parameter n is the radial location with n - 0 
being the center of the pipe. The traces show a laminar-like flow during 
acceleration and a turbulent-like flow during deceleration. With increased 
Re max (5830, Fig. 3-7b), the turbulence persists into the acceleration 
portion of the cycle. From these traces, one would expect transition to 
extend over a broad range of Remax* This supposition is consistent with 
measurements taken by Ohmi et al. (1982) where a wide range in Re m ax was 
observed between laminar-like and turbulent-like oscillating flows. The 
signal in the first part of the cycle shows slightly stronger fluctuations 
than that in the second part, probably because of flow separation in the 
asymmetric apparatus. At the same dimensionless frequency (Re^ _ 7.30) 
but a higher amplitude Reynolds number (Re max ■= 5830) the turbulent 
fluctuations are stronger than at Remax “ 2070 (Figure 3~7(a); they start 
prior to deceleration during the high velocity part of the cycle. 
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The ensemble-averaged velocity is lower than the laminar velocity at 
the center of the pipe and higher than the laminar velocity near the pipe 
wall. This indicates that the profile becomes flatter, which is 
characteristic of turbulent velocity profiles. This observation is 
consistent with the work of Ohmi et al. (1982) who found that the velocity 
profile during the laminar part of the cycle agrees well with the 
theoretical laminar solution (see Section 3.1) and that the velocity profile 
during the turbulent part agrees well with the 1/7 power law for steady 
turbulent pipe flow (see e.g. Schlichting 1979, p. 597-602). 

Experimental observations of transition . Figure 3-8 shows observations 
of transition. Ohmi et al. (1982) studied forced oscillations of a gas in a 
straight pipe. They did not state their criterion for transition. Iguchi 
et al. (1982) observed free oscillations of a liquid column in a U-tube. 

They chose the laminar-to-transitional flow line to be where the amplitude 
of oscillation began to deviate from that of the laminar solution discussed 
earlier. The transitional- to- turbulent flow line is where measured 
amplitudes of oscillation began to agree with predictions based on 
1/7th-power turbulent pipe flow profiles. Park and Baird (1970) observed 
the decay of free oscillations of liquid in a manometer. They calculated 
cycle maximum wall shear stress from the observed amplitudes of oscillation 
using two methods. The first was based on a laminar flow prediction, the 
second on the turbulent 1 /7th-power profile. They assumed transition to 
occur at the amplitude where the maximum wall shear stress, calculated from 
the laminar prediction, exceeded that calculated from the 1/7th-power 
profile. They attribute their data scatter to column end effects in the 
liquid column which are a function of L 2 u>/v where L is the length of the 
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liquid column. These end effects may be important in liquid piston Stirling 
engines. Sergeev (1966) did not state his criterion, but. since he used 
aluminum particle flow visualization, we surmise that transition was 
observed through the transparent walls as a change in the flow structure. 

His work was in straight-tube, forced-oscillation flow. Grassmann and Tuma 
(1979) used an electrolytic technique to observe turbulent fluctuations in 
wall mass transfer rate, thus locating transitional behavior. 

Merkli and Thomann (1975) studied the transition from laminar to 
turbulent oscillating flow at dimensionless frequencies (Re^ - 
2500. ..4000) beyond the range presented in Figure 3-8. Their apparatus was 
a piston oscillating in a resonance tube. They noticed a secondary weak 
vortex motion outside the oscillating boundary layer. The tube was operated 
near resonance. The secondary vortex observed by Merkli and Thomann may be 
related to the vortex street that Sobey (1985) observed in steady and 
oscillating channel flow. Sobey also predicted the vortex street 
numerically and considered it a result of shear-layer instability. 

Dijkstra (1984) observed, in the transitional regime, that water which 
was initially in the tube would remain laminar, while water entering the 
tube would flow in a turbulent pattern and would remain turbulent in the 

tube. 

Though the transition predictions differ with criteria, the researchers 
agree that transition Re max increases with Re u . The sequence of the 
transition predictions is consistent with the criteria used. Iguchi’s et 
al. (1982) lower line is based on the first sign of deviation from laminar 
behavior. The Grassmann and Tuma (1979) criterion is based on fluctuations 
at the pipe surface. These fluctuations are likely to occur at about the 
same Re M as that established by Sergeev (1966) based on the onset of 
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turbulent motion of particles. Finally, the Iguchi et al. (1982) upper line 
is based on agreement with the predictions from the 1/7th-power law. It 
agrees with the observations by Ohmi (1982) in forced oscillation of gas in 
a straight pipe. Equations for these observations of transition are listed 
in Appendix D. 

Theoretical prediction of transition . Von Kerczek and Davis (1972) 
used the energy method to predict a lower bound for the instability of 
Stokes layers on a flat plate. Figure 3-9 shows their results in the form 
of the critical value of the Reynolds number Re max below which the 
oscillating flow cannot go unstable. This lower bound underpredicts the 
transition Reynolds number by roughly one order of magnitude. The trends of 
the two curves, however, agree well. 

Cayzac et al. (1985) presented predictions for the lower bound of 
stability in oscillating pipe flow. Like Von Kerczek and Davis, they found 
good qualitative agreement in the trends but a quantitative discrepancy of 
one order of magnitude. For steady pipe flows, they predicted the lower 
bound for instabilities to be Rep » 750 while Rep ■ 2000 is the well-known 
experimental value. 

The theoretical prediction of the lower bound of instability by the 
energy method apparently does not yield results that are practically useful. 
For other theoretical approaches to the stability of oscillating flows see 
Davis (1976). 
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Figure 3-9: Theoretical prediction of transition 

and experimental observations 
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3.5 Turbulent Flow 

Though much remains to be learned about the effect of oscillation on 
turbulent flows, significant findings have been documented. 

Quasi-steady approach . Dijkstra (198M) and Vasiliev and Kvon (1971) 
applied turbulence models taken from steady flow to oscillating and 
pulsating pipe flow, respectively. Kirmse (1979) compared Vasiliev and 
Kvon's model to his experimental data and concluded that their quasi-steady 
prediction was not adequate. 

The concept of an oscillation-sensitive and a quasi-steady turbulent 
regime was introduced by Ohmi et al. (1982, p. 536). In their experiment, 
the crank radius of the slider-crank mechanism driving the flow was 
increased thereby increasing Re max at constant Re u . This increased the 

crank radius to connecting rod length ratio (r/l). Figure 3—1 0 shows plots 

« •• 

of the displacement (x), the velocity (x) and the acceleration (x) of a 
slider-crank mechanism. For r/i, - 1/3 the acceleration, x, shows large 
deviations from a sine function though the displacement, x, is nearly 
sinusoidal. The turbulence structure is known to be sensitive to 
acceleration and deceleration. It is, therefore, not possible to decide 
whether deviations from sinusoidal, quasi-steady flow are due to changing 
the drive geometry or due to the increase in the ratio of Re max / Re^. 

Ohmi et al. (1982) claimed the latter when they determined the line Re max « 
2800 /Re u above which they predict turbulent flow to be quasi-steady. 

They applied a quasi-steady 1 /7-power law velocity profile and observed that 
it agreed well with their measured velocity profiled above this line. 
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Fluctuating eddy viscosity model . Kita et al. (1980) found that the 
Reynolds stresses In pulsating flow could not be adequately described using 
a constant eddy viscosity. Therefore, they proposed a fluctuating eddy 
viscosity model. Their five-layer model shows good agreement with 
experimental results for pulsating flows; it wa3 not tested for oscillating 

flows. 

Empirical pressure drop correlations . Taylor and Aghili (198*0 took 
data of pressure drop in oscillating flow through pipes of finite length 
They plotted their data as a function of time-averaged Reynolds number for 
various values of A/d. Their "friction factor" is based on a time-averaged 
pressure drop and a time-averaged absolute value of the mean velocity. The 
data indicate an increase of the friction factor by a factor of four over 
the steady, unidirectional case. 

From the previous discussion of pressure drop in laminar oscillating 
pipe flow, it is clear that the acceleration and deceleration of the working 
fluid influence the pressure drop in oscillating flow. Only if the results 
were correlated in terms of a dimensionless frequency such as Re^, could 
they be applied appropriately. Taylor and Aghili* s (1984) paper does not 
give sufficient data to calculate Re u . 

3.6 Entrance and Exit Losses 

Thus far, the discussion of pipe flow was restricted to fully developed 
flow or, equivalently, infinitely long pipes. For finite pipes, entrance 
effects may be important. The viscous losses are higher in the entrance 
length of a pipe. Because of the spatial acceleration of the fluid, 

u*Vu term in the momentum equation is not zero. 
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Goldberg (1958) showed that the hydrodynamic and the thermal entrance 
length in hydraulically and thermally developing steady laminar flow is a 
function of Reynolds number. Based on this observation, Charreyron (1984) 
suggested that the entrance length varies over the cycle. Peacock and 
Stairmand (1983) hypothesized that the entrance length in laminar 
oscillating flow will be shorter than in unidirectional, steady flow. Since 
the velocity profiles tend to be flatter in oscillating flow, the velocity 
profile of the oncoming flow may not have to change much in the entrance 
region. Their hypothesis so far was not supported by experiments. 

Disselhorst and van Wijngaarden (1980) studied separation near the 
entrance of a tube under conditions of acoustic resonance. They found that 
separation did not occur for high Strouhal numbers. During the inflow, a 
boundary layer forms on the pipe walls; during the outflow, a jet emerges 
from the pipe. They observed and predicted that vortices that were shed 
during the outflow would interfere with the formation of the boundary layer 
during inflow. 

3.7 Compressibility Effects 

Compressibility effects in Stirling engine heat exchangers may be due 
to 

(1) pressure recovery when the working fluid exits from a duct or a 
heat exchanger tube into a large volume such as the cylinders, 

(2) high Mach numbers in heat exchangers, 

(3) travelling shock waves caused by the interaction of pressure waves 
generated by the motion of the pistons, 

(4) finite speed of pressure propagation. 
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Pressure recovery . The effect of the increase in pressure on the gas 
density can be approximated by an adiabatic expansion. This is commonly 
done by a compressibility correction factor ("expansion factor") for flow 
through nozzles or sudden expansions (e.g.. Perry et al. 1973. p. 5-11). 

High Mach numbers . If the Mach number of the fluid exiting from a heat 
exchanger is one, a shock will form at the exit and the flow will be choked, 
thereby limiting the mass flow rate. It will be shown later that these high 
Mach numbers are not expected for Stirling engines discussed in this report. 

Travelling shocks . Research on non-linear acoustics in a resonance 
tube shows that travelling shocks can form near resonance conditions at Mach 
numbers less than one (M max » 0.2) (Merkli and Thomann 1975 a and b). This 
was predicted theoretically by Jeminez (1973). Resonance conditions, 
however, are unlikely in Stirling engines because the flow length is 
typically much less than half the acoustic wavelength (see Appendix E). 
Travelling shocks can develop without resonance, however, by the following 
mechanism. 

When a piston moves towards top dead center, it causes a compression 
wave to travel from the piston face. This wave travels through the 
cylinder, ducts and heat exchangers. The high-pressure part of a wave 
travels faster than the low-pressure part of the wave. Therefore, the 
pressure wave steepens and may form a shock (Shapiro 1 95-4 , p. 9^9-950). 

This process is affected by friction and heat transfer only through the 
change in the local velocity of sound due to change in local fluid 
temperature. A simple analysis of this type of shock formation is carried 
out in Appendix E. This worst case estimate assumes that the working fluid 
motion is driven by a piston with a velocity amplitude of M max - 0.15 (the 
maximum value found in the data base, see Section 7). The fluid temperature 
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is assumed to vary linearly from T c - 350K to - 1050K throughout the 
engine and the frequency is the highest one found in the data base (SPDE-D), 
f - 105 Hz. Appendix E explains why this may be considered a worst case. 

The results of this worst case show that formation of shocks in the 
Stirling engines in this data base is not expected. The analysis is limited 
in that it only includes right-travelling waves and does not take into 
account the interference of right- and left-travelling pressure waves from 
the opposing pistons of an engine, which would distort the Mach lines used 
to predict shock incipience. 

A more complete analysis of compressible flow in Stirling engines was 
carried out by Organ (1982b) and modifications were proposed by Taylor and 
Aghili (1984). Organ's analysis predicts that there is no incipient shock 
for the special case of the air-charged engine discussed in his paper. This 
section shows that incipient shock formation is unlikely in the engines in 
the data base. 

Finite speed of pressure propagation . Pressure changes propagate with 
finite speed, i.e., the sum of the velocity of the piston causing the 
pressure change and the velocity of sound. Therefore, a phase lag is 
expected between the pressures at both ends of a heat exchanger. If this 
phase lag is a small part of the cycle, pressure propagation throughout the 
heat exchangers may be considered instantaneous. 
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4. HEAT TRANSFER IN OSCILLATING PIPE FLOW 


Since convective heat transfer depends on the velocity distribution, 
the similarity parameters for fluid mechanics affect the heat transfer: 
Re raax , Re <u» a R or A/d * In additlon > similarity parameters from the energy 
equation are important: the Prandtl number, Pr, and for high speed flows, 

the Eckert number, Ec.the Mach number, M, and fyj/d as a measure of the 
thermal entrance length. The wall temperature, T w , the bulk fluid 
temperature at the hot end, T^, and at the cold end, T c , also control heat 
transfer. 
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4.1 Qualitative Considerations 


Limited fluid displacement . The limited streamwise displacement of 
fluid in Stirling engines may limit heat transfer from the heat exchanger 
wall or matrix as discussed in Section 2.1. Figure 4-1 shows fluid 
trajectories based on the plug flow assumption. While 8,^ is the 
characteristic streamwise length for the heat transfer problem, the flow 
length Jl is used in this study because i. » in most Stirling engines and 
in many cases 8,^ n °t documented. 



Figure 4-1 : Fluid displacement in heat exchangers. 
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Under the plug flow assumption, the Ar values indicate: 

Ar < 1 Some fluid does not leave the heat 

exchanger. 

Ar - 1 All fluid moves into and out of the 

heat exchanger. 

Ar > 1 Some fluid moves through, some fluid 

moves in and out of the heat exchanger. 

Assessing plug flow, Organ (1975, p. 1016) argued that, in the absence 
of axial mixing, the fluid that never leaves the heat exchanger in the case 
of Ar < 1 could be considered additional dead volume. In turbulent flow, 
the presence of axial mixing may invalidate Organ’s model. 

The velocity profiles for laminar flow in Figure 3~2 and in Appendix C 
show that plug flow is not a good assumption unless Re u is very high. In 
turbulent flow, plug flow may be a good assumption but axial mixing will 
become important. The velocity profiles show that, in fully developed 
oscillating pipe flow, the fluid displacement varies with radial distance. 
Therefore, Organ's (1975) concept of moving temperature wavefronts is only a 
rough approximation for turbulent flows. In fact, the radial variation in 
fluid displacement leads to enhanced axial heat transfer as discussed in 
Section 4.2. 

Two temperature driving potentials . Typical boundary conditions for a 
tubular heat exchanger in a Stirling engine are sketched in Figure 4-2. 
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Th > T w > T c 


Figure 4-2: Boundary conditions for heat transfer in a heater. 

During the hot blow, gas from the expansion space (at temperature T^) will 
stream through the tube into the regenerator. The laminar velocity profiles 
in Figure 3-2 show that, at higher values of Re^, the fluid near the tube 
wall may flow in one direction while the core moves into the opposite 
direction (Re^ - 30, crank angles 150° and 180°). For those parts of a 
cycle where counterflow occurs, an unusual heat transfer situation exists: 
the wall gives off heat to the fluid near the wall which enters the tube at 
T c from the regenerator. The fluid coming from the expansion space receives 
heat from the near-wall-fluid in a counter-flow heat exchange process. 

In this case, the heat transfer has two driving potentials, (T w - T c ) 
and (T h - T c ) where (T w - T c ) > (T h - T c ) . It appears that backflow near 

the wall may have a significant effect on heat transfer. Its influence 
would depend on the degree of mixing within the near wall shear layer. The 
concept of two driving potentials for heat transfer in oscillating flow is 
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also important because of the augmentation of heat transfer in the axial 
direction discussed in the next section. 

A situation of two temperature driving potentials also exists in film 
cooling of gas turbine blades, where the temperature differences between the 
hot gas in the free stream and the injected cooling air determines the heat 
transfer in addition to the temperature difference between the free stream 
and the blade surface. A concept like the film cooling effectiveness may 
prove useful for the interpretation of experimental results on heat transfer 
in oscillating flow. 
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4.2 Axial Heat Transfer In Laminar Oscillating Flow 

This section discusses how energy transport due to an axial temperature 
gradient is enhanced by flow oscillation. Note that this augmentation does 
not depend on turbulent cross-stream mixing and that the analysis only holds 
for laminar pipe flow. 

Physical interpretation . When fluid oscillates in a duct in the 
presence of an axial temperature gradient, large oscillating temperature 
gradients normal to the flow direction are generated. During half of the 
cycle, cold fluid at the core of the duct extracts thermal energy from the 
hot fluid in the boundary layer. During the other half of the cycle, hot 
core fluid heats up cold fluid near the pipe wall. Thus heat is transferred 
from the hot end to the cold end of the duct. For further physical 
interpretation, see Kurzweg (1985b, p. 298). 

Note that this physical interpretation , as well as the analyses 
discussed below, hold for laminar flow. In turbulent flow, convective 
eddy-transport in the cross-stream direction may reduce the oscillating 
cross-stream temperature gradients that are necessary for the augmentation 
of axial heat transfer. 

Analysis . Watson (1983) analyzed diffusion in laminar, oscillating, 
fully developed flow in impermeable ducts of arbitrary, uniform 
cross-section. The analogous heat transfer situation is sketched in Figure 
4-2 except that the walls are adiabatic. The pipe connects a high 
temperature reservoir (Tp) and a low temperature reservoir (Tl). Assume 
that end effects are negligible and that Ap < 1 . In a stagnant fluid, 
assuming that there is no convection, the axial heat flow rate is due to 
conduction : 
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Watson (1983) showed that the flux through a duct in which an Incompressible 
fluid oscillates in laminar flow can be expressed similarly as 



where the apparent thermal conductivity 
k eff - (1 + Au) 

depends on the augmentation coefficient Au, which is a function of Re^, Pr 
and the dimensionless tidal volume 



The augmentation coefficient is 

Au - Au (Re u , Pr, V T ) - V T 2 ? (Re^, Pr) 
is given in Appendix F for values of Pr and Re u that are relevant for 
Stirling engine heat exchangers. Kurzweg (1985a, p. 461) provides a plot of 
Au vs. /Re^ for Pr - 0.1, 1., 10. He develops an asymptotic solution for 
Re^ Pr < tt. Watson's (1983) exact solution wa3 confirmed experimentally 
by Joshi et al. (1983) for the case of a circular pipe. Gedeon (1986) and 
Kurzweg (1985b) treat the case of flow between parallel plates. Kurzweg 
(1985b) calculates the augmentation of axial heat transfer in the case where 
the walls are diabatic, i.e. they undergo transient temperature fluctuations. 
The walls contribute to heat transfer by absorbing thermal energy from the 
hot fluid and giving off thermal energy to the cold fluid. Kurzweg* s 
(1985b) analysis shows that the augmentation of axial transport is a 
function of Pr, Re^, V T , and the ratios of fluid to wall conductivity and 
thermal diffusivity. 
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4.3 Experimental Data 


Hwang and Dybbs (1980 and 1983) presented experimental heat transfer 
results for oscillating flow in a tube. Figure 4-3 shows a schematic of 
their experimental facility. Air was moved through a heater and cooler by 
an oscillating piston into an adiabatic cold space. The piston moved 
approximately sinusoidally. As it moved to the right, the air expanded in 
the cooler and into the heater. The heat transfer was calculated by 
measuring the temperature rise of the cooling water on the outside of the 
test section. Thermocouples measured the (apparently time averaged) gas 
temperatures at both ends of the cooler and the wall temperatures in the 
cooler pipe. 



cold space cooler heater cylinder piston with 

(test section) slider crank 

mechanism 


Figure 4-3: Schematic of Hwang and Dybbs' (1980) experiment 

A Nusselt number was then formed by the use of the log-mean-temperature 
difference, the measured average heat flux, and the pipe diameter. The 
results are replotted in Figure 4-4 by the present authors. 
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Figure 4-4: Hwang and Dybbs* (1980) data 
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Hwang and Dybbs (1980 and 1983) plotted Nu vs. Re oa (Re 03 - Re max ^ 
(2it)) with A m - 2/Ar as a parameter. In this experiment only two of the 
parameters Ar, Re max and Re^ could be specified independently. Ar is 
physically not significant since for most data points the gas will be blown 
completely through the heat exchanger (Ar > 2) , thus Re max and Re^ remain 
as the only significant parameters. 

The data in Figure 4—4 show certain trends: 

(a) For most values of Re w , the data show a small slope for low 
values of Re max , then the slope increases with increasing 
Re^. 

(b) The data points with the higher slope tend to lie in one 
common band. 

(c) The individual data points typically scatter within a band 
that is wider than two standard deviations of a single 

data point; i.e., the scatter of the data cannot be accounted 
for by the random error of single data points. 

The lower Re max data points of each curve may represent predominantly 
laminar flow patterns over a cycle while the parts of the curves with 
greater slope may represent turbulent flow during parts of the cycle. This 
is argued by Hwang and Dybbs analogous to unidirectional steady flow where 
Nu a Re^3 in laminar flow and Nu a Re®*® in turbulent flow. In view of 
statement (c) on the scatter of data no further conclusions are drawn. 

The heat transfer in Hwang and Dybbs' experiment consists of two 
different processes: 

(a) blow of hot air through the cooler pipe into the cold space. 
During this blow the heat transfer may be approximated by 

q" « h A_ 

x lm 

where ATi m is the log mean temperature difference 

and h - — kg - conductivity of the gas 

This can be computed if the instantaneous heat flux and 
log-mean temperature difference (LMTD), or at least the time 
average values for this blow is known. 
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(b) blow of cool air through the cooler pipe into the heater. 

Again, a value of Nu' based on the time-average heat flux and 
LMTD for this part of the cycle could be calculated. 

In Hwang and Dybbs' experiment, however, only the heat flux averaged 
over the whole cycle and the LMTD based on temperatures averaged over the 
whole cycle are known. Thus the Nusselt number calculated is neither 
representative of the hot blow nor of the cold blow, but is an average of 
the two. The heat flux that Hwang and Dybbs used to calculate a Nusselt 
number may be due to heat transfer during the hot blow only. The LMTD which 
is apparently based on time-averaged temperatures probably underestimates 
the temperature differences during hot blow. Thus it is difficult to draw 
further conclusions from this data set. 

Apparently, none of the transition predictions fit Hwang and Dybbs' 
data (Figure 4-5). Possible explanations for this discrepancy are: 

(1) The tube was very short (4/d - 28); therefore the transition 
criteria which were presumably developed for fully developed 
oscillating flow may not apply. 

(2) The first data point on each curve may be taken from relative 
amplitudes that are too small to measure the heat transfer 
coefficient accurately (Ar - 1.8). 

Neither criterion may be sufficient to explain the change in 

time-averaged heat transfer for the Hwang and Dybbs experiment because the 
heat transfer is a function of: 

(a) the location where turbulence occurs (e.g., near the wall or 
in the core). 

(b) the relative time during which laminar and turbulent flow are 
present during the same cycle. 
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Figure 4-5: Comparison of prediction and observation of transition. 


Iwabuchi and Kanzaka (1982) studied heat transfer in oscillating flow 
in a test facility that was designed to obtain results for the design of a 
specific prototype engine. A rough estimate indicates that they operated in 
the laminar regime (11 < Re^ < 83, 145 < Re max ^ 733). They did not 
correlate their results in terms of a Reynolds number and a dimensionless 
frequency but they studied the dependence of the heat transfer on parameters 
such as rpm, mean pressure and phase difference between the two opposing 
pistons. Their heat-transfer results are, therefore, not generally 
applicable. They made the observation, that the choice of the phase 
difference (90° or 180°) did not change the heat transfer provided that the 
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Schmidt analysis was used to calculate the mass flow of gas through the heat 
exchanger . 
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5. FLUID MECHANICS AND HEAT TRANSFER IN REGENERATORS 


In the complex geometries of regenerators, flow separation is expected. 
Therefore, analytical methods are not applicable and experimental results 
are required. 

5.1 Steady Flow 

Stacked, woven wire screens, randomly stacked metal fibres, folded 
sheet metal, metal sponge and sintered metals are used as regenerator 
matrices. Steady flow through these matrices has been studied extensively 
in porous media research and in studies related to gas turbine and Stirling 
engine regenerators. 

Flow Regimes in Porous Media Flow 

Dybbs et al. (1984) classified the steady flow regimes in porous media 
according to a Reynolds number based on the average pore diameter. 

Re < 1 
1 < Re < 10 
10 < Re < 175 
175 < Re < 250 
250 < Re < 300 
Re > 300 

The Reynolds numbers separating these regimes are expected to hold for only 
small values of Re w and only for materials similar to those used by Dybbs 
et al. (1984), but they represent the only known work on transition that 


Darcy flow regime 

Boundary layers begin to develop on the pore wall 
laminar flow 

separated laminar flow, vortex shedding 
separated flow with random wakes 
turbulent flow 
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would apply to regenerator matrices. The complex and irregular geometries 
of porous media and the complex flow patterns in flows with Re > 175 make it 
difficult to predict pressure drop and heat transfer with numerical methods. 
Therefore experimental results are required for the prediction. 

The Permeability Model 

Beavers and Sparrow (1969) showed experimentally that the pressure drop 
for steady flow of an incompressible fluid through a porous material can be 
described by 

-dp/dx - apV + bpV 2 
where: 

a, b are constants to be determined; 

V is the superficial velocity; 



A is the frontal area of the porous material. In this equation, the 
first term accounts for purely viscous pressure drop as it is described by 
Darcy’s law: 

-dp/dx - pV/k 

where the permeability, k, is a constant for each porous material. 

Beavers and Sparrow were able to reduce pressure drop data for foamed nickel 
(FOAMETAL by General Electric) of various geometries and porosities and for 
wire screens by using the square-root of the permeability as the 
characteristic length. The Reynolds number and friction factor were defined 



r f_ lEl ^L, 

^ dx' W 

Then the dimensionless pressure drop equation is: 
f - 1/Re + C 


The values for C were found to be close to C « 0.074 for FOAMETAL and wire 
screen, while a specimen of FELTMETAL (randomly stacked fibers, Huyck Mfg. 
Co.) had a higher value of C - 0.132. There are free fiber ends in 
FELTMETAL which are not found in the other materials. This structural 
difference may be the reason for the difference in C. 

To see whether the results obtained by Beavers and Sparrow are relevant 
for Stirling engine flow conditions, the present authors calculated the 
Reynolds number for the wire screen specimen, which Beavers and Sparrow used. 
This Reynolds number, comparable to Re max , ranged from about 64 to 10* 4 . 

In order to use the friction factor correlation 
f - 1/Re + 0.074 , 

the permeability must be known. This quantity is obtained experimentally. 
The pressure drop of creeping flow is measured and the permeability is 
calculated from the experimental data using Darcy’s law as given above. 

Dybbs et al. (1984) studied steady flow in porous media. They observed 
that at higher velocities than those characteristic of Darcy flow, the 
friction factor is a function of a Reynolds number based on the hydraulic 
diameter and on the ratio of length to diameter of a typical pore. Their 
explanation was that the flow must develop anew in every pore. For the 
prediction of pressure drop and heat transfer, this may mean that the 
length-to-diameter ratio of a typical pore should be included as an 
additional geometric parameter. The permeability cannot represent the 
effect of the pore shapes properly because it is measured in the 
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non-inertial Darcy regime. This is one reason why only geometrically 
similar materials have the same C-value as observed by Beavers and Sparrow 
(1969). 

Macdonald et al. (1979) reviewed research on flow in porous media 
consisting of spherical beads, cylindrical fibers, sand, gravel, and many 
others. They obtained an equation that predicts pressure drop for a wide 
range of materials within ±50$. 

Joseph et al. (1982) showed that the pressure drop in a bed of packed 
spheres can be related to the drag on a single sphere. They also lent 
theoretical support to the velocity-squared term proposed by Beavers and 
Sparrow (1969). 

Beavers et al. (1973) studied the influence of the shroud bounding 
porous media. They found that the pressure drop in beds of packed spheres 
was influenced by the walls for shroud-diameters as small as MO sphere 
diameters. 

Friction Factor Correlations 

Kays and London (1984, p. 149) provide widely used correlations for 
flow through stacked screens. Walker and Vasishta (1971) present 
experimental data for dense-mesh wire screens. Miyabe et al. (1982) provide 
additional experimental data for flow through stacked screens, Takahashi et 
al. (1984) present data for foamed-metal matrices. Chen and Griffin (1983) 
derived an empirical equation based on the Kays and London (1964) data that 
more accurately fits the data. A comprehensive review of regenerator 
pressure drop and heat transfer correlations for steady flow was prepared by 
Finegold and Sterrett (1978). The disadvantage of these correlations is 
that they do not collapse data as well as the correlation proposed by 
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Beavers and Sparrow (1969) because porosity remains as a parameter in 
addition to the Reynolds number. The advantage is that they were obtained 
for gas flows so that some compressibility effects such as reduced effective 
flow area, which is discussed below, are Included. 

Compressibility effects in flow through screens . Organ (1984) 
recommended that the Mach number be included as a correlating parameter for 
the prediction of pressure drop in regenerators. His suggestion is based on 
earlier work by Benson and Baruah (1965) and Pinker and Herbert (1967) where 
the flow of air through a single screen was investigated. 

Benson and Baruah (1965) stated in their conclusions (p. 458) that "The 
resistance coefficient of a gauze is a function of Mach number, Reynolds 
number and solidity," as quoted by Organ, solidity - 1 -porosity. In the 
same paper they plot the resistance coefficient of screens (with porosities 
of 0.47... 0.67) versus downstream Mach numbers. For low Mach numbers (M < 
0.1) Benson and Baruah's (1965) results show that the resistance coefficient 
is virtually independent of Mach number. Benson and Baruah (1965) also 
state in their conclusions (p. 458) that a "gauze behaves like a crude 
nozzle. There is a considerable entropy change in a flow across it." 
According to Pinker and Herbert (1967, p. 16 and p. 7) this "nozzle" 
develops at Mach numbers of the oncoming flow of M - 0.2 in the duct and of 
M - 0.38 for a porosity of 0.53. The duct Mach number required for choking 
increases with porosity. 

Effective Flow Area in Compressible Flow 

To calculate the mass flow rate through a nozzle, the flow area is 
multiplied by an expansion factor (see e.g.. Perry et al. 1973, p. 5-11). 

It accounts for the reduction in effective flow area due to the adiabatic 
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expansion from the upstream to the downstream pressure. It Is expected that 
a similar effect can be observed when flow passes through a single screen. 

It is not likely that this effect will be significant for flow through 
stacked screens because the pressure differences across individual screens 
will be small and heat transfer from the matrix to the gas will be possible. 

Choking and Pressure Drop in Compressible Flow 

Beavers and Sparrow (1971) calculated the conditions under which the 
flow into a porous material would be choked based on isentropic flow. 
According to Benson and Baruah (1965, p. ^58), a considerable change in 
entropy takes place when a gas passes through at higher Mach numbers. This 
entropy change may not be significant for the low Mach numbers encountered 
in Stirling engines, however. 

Beavers and Sparrow (1971, p. 1856) also analyzed gas flow within a 
porous matrix using /iT (the square-root of the permeability) as the length 
scale. They calculated the maximum length J^ax of porous material before 
the exit flow would be choked. A regenerator operated at high Mach number, 
low Reynolds number and with a high iSk is most likely to choke. Even in 
the extreme case of the cold blow through the ML23 regenerator H. m ax /y/ ^ > 
lOfc/i/Jr. Therefore choking is not expected In regenerators. 

This analysis is tentative because: 

(1) The permeabilities were calculated from the comparison of 
permeability and hydraulic diameter for one specimen that Beavers and 
Sparrow (1969) used; to calculate the permeabilities for other wire-mesh 
regenerators a linear relationship between hydraulic diameter and the 
square-root of the permeability was assumed. 
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(2) the graphs used to calculate l ma x / ^’ ln Beavers and Sparrow (1971) 
were plotted. . . 

(a) for air only. 

(b) for discrete values of Reynolds number. 

(c) for discrete values of porosity. 

Heat transfer correlations . Kays and London (1984, p. 149), Walker and 
Vasishta (1971), Finegold and Sterrett (1978), and Miyabe et al. (1982) 
provide heat transfer correlations for packed woven wire screens. These 
correlations are derived from experimental data that were obtained by the 
single blow technique described in Kays and London ( 1 98 1 * , pp. 15*1-155). 
Takahachi et al. (1984) calculated heat-transfer correlations from 
oscillating flow data. The use of A” as the characteristic length instead 
of the hydraulic diameter may eliminate the porosity as the independent 
parameter in heat transfer correlations. This remains to be shown 
experimentally . 

5.2 Unsteady Flow and Heat Transfer in Regenerators 

Galitseiskii and Ushakov (1981) studied heat transfer augmentation due 
to reversing pulsatile flow in porous media. They observed a resonant 
phenomenon when the diameters of the stationary vortices (generated by the 
flow over the elements of the porous structure - e.g. the wires of the 
screen) are comparable to the dimensions of the secondary vortices generated 
by the flow oscillations. The experiments were performed for 0.6 < Re max < 

6 and 6.5 < Re u < 100, a domain that does not seem applicable to 
regenerator flow conditions. 

Kim (1970) studied oscillating flow through a regenerator consisting of 
packed spheres. The reduction of his experimental data was based on the 
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assumption that the momentum equation could be considered quasi-steady i.e., 
the pressure drop is in phase with the mass flow rate. He showed by 
order-of-magnitude estimates that this was justified (p. 130). He 
calculated friction factors for gas flow through a bed of packed spheres 
that were approximately 20$ higher than the corresponding steady values. 

The range of frequencies covered less than a decade. Since the stroke of 
his apparatus was fixed, the dimensionless frequency of his experiments was 
proportional to the peak Reynolds number. Therefore, his data do not allow 
conclusions about the frequency dependence of the friction factor. His 
results do, however, suggest that fluid friction increases due to flow 
oscillation. Because the friction factor is higher in oscillating flows 
than in steady flow, it must depend on frequency (Rem). Whether the same 
is true for wire-screen regenerators must be answered by experiment. 

5.3 Regenerator Theory 

Stirling engine simulation typically uses a discretization of the 
engine volume. If the regenerator is subdivided in several control volumes 
it is not necessary to use regenerator effectiveness to predict heat 
transfer to and from the regenerator (Walker 1980, p. 140-149). Given an 
adequate regenerator theory, however, measured effectiveness values can be 
used to determine volume-averaged heat transfer coefficients of regenerators. 
This method has been used successfully for gas-turbine regenerators (Kays 
and London 1984, p. 154-155). 

Rice et al. (1985) proposed a theory that takes into account the 
relatively short blow periods in Stirling engines. Willmot and Hinchcliffe 
(1976) and Harness and Neumann (1979) proposed methods to take "holdup" into 
account, i.e., the fact that some fluid may not leave the regenerator by 
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convection (A R < 1). These approaches seem more promising than the 
application of standard regenerator theories (e.g. Shah 1981, p. 721-763 or 
Kays and London 198H, p. 79-89) to Stirling engine conditions, e.g. by Rice 
et al. (1983) and Miyabe et al. (1982). 
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6. STIRLING ENGINE DATA BASE 


Representative values of the similarity parameters describing the 
operating conditions of the heat exchangers of 11 Stirling engines are 
obtained to form a Stirling engine heat exchanger data base. 


6.1 Evaluation of Similarity Parameters 


Schmidt analysis . The instantaneous mean mass velocity pu m is 
calculated using the Schmidt isothermal analysis (Uriel! and Berchowitz 
1984, Appendix A1 , p. 152). The following assumptions are made for the 


analysis : 

1. The engine consists of three isothermal volumes: 

The hot space (at T^) consists of expansion space, connecting 
duct, heater, 1/2 heater-regenerator duct. The regenerator 
volume (at T r ) consists of 1/2 heater-regenerator duct, the 
regenerator and 1/2 regenerator-cooler duct. The cold space 
(at T c ) consists of compression space, connecting duct, 
cooler, 1/2 cooler-regenerator duct. T r is the mean effective 
regenerator temperature (Urieli and Berchowitz 1984, Appendix 
A4, p. 158) 


Tr 


T h - T c 
ln(T h /T c ) 


2. The pressure is uniform throughout the engine. 


3. The ideal gas law holds. 

4. The variations of the compression and the expansion volume are 
sinusoidal . 


5. There is no working gas leakage. 

6. The engine operates at constant speed in a cyclic steady 
state. 

7. The kinetic and potential energies of the gas are negligible. 


Isothermal analysis . The isothermal analysis is based on the same 
assumptions as the Schmidt analysis except point 4. The variations of 


compression and expansion volumes are calculated from the kinematics of the 
drive mechanism. Appendix A shows Reynolds number variations calculated for 
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a rhombic drive engine (GM GPU-3) with Schmidt (Figure A-3) and isothermal 
analysis (Figure A-1) (Youssef 1986). The two analyses differ in the shape 
of the Reynolds number transients. The difference in amplitude of 
approximately 20 percent is acceptable because this data base is only used 
to establish approximate ranges of heat exchanger operating conditions. 

Reynolds and Mach number . The thermodynamic analysis is used to 
calculate the changes of mass in various control volumes of the engine. 

From the mass changes in the control volumes, mass flow rates are 
calculated, and from these Reynolds and Mach numbers. u m , the mean 
velocity, is the velocity averaged over the flow area, e.g. over the 
cross-sectional area of a duct. Figure A-1 shows the variation of the 
Reynolds numbers at the inlet and outlet of a GPU— 3 heater calculated by the 
isothermal analysis. The values at inlet and outlet differ because working 
gas is stored in the heat exchanger during the compression and released 
during the expansion part of the cycle. This changes the mass flux term pu m 
in the Reynolds number. Therefore, the maximum Reynolds number at each 
operating point is represented by a range of values. This range represents 
the variation throughout the heat exchanger due to asymmetric volume 
variations in storage and release of working gas during compression and 
expansion (see Appendix A). The dynamic viscosity is evaluated at the 
temperature of the corresponding control volume, the density is calculated 
from the ideal gas law, and the speed of sound is evaluated assuming 
constant specific heats, a - /WT. The hydraulic diameter used to calculate 
Re max for non-tubular geometries (rectangular ducts) is defined as the ratio 
of volume to heat transfer area of the heat exchanger. 

Dimensionless frequency Re,.' and relative amplitude Ar . The 
dimensionless frequency is evaluated at the mean pressure and the 
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temperature of the respective heat exchanger. For regenerators, the 
hydraulic diameter is used as a length scale. The relative amplitude of 
fluid motion, Ar, is evaluated by integrating the mass flow at the heat 
exchanger inlet or outlet between two flow reversals. For the evaluation of 
the displacement of a fluid element between flow reversals the fluid is 
assumed to move into and out of the heat exchangers as a plug. 

6.2 Documented Engines 

The choice of documented flow conditions was restricted by the data 
available. Table 6-1 lists (alphabetically by abbreviation) the name of the 
engines included in the data base, a brief description of the operating 
conditions and the reference used. The amplitudes of the pistons and 
displacers of the free pi3ton engines were not obtained from the analysis 
but experimental values were taken from the references. 

The SPDE "design" point (SPDE-D) was calculated from phase angle and 
amplitudes, measured in the trial run represented by "SPDE-T" but using 
design pressure and frequency. A more realistic approximation to the 
intended operating conditions was calculated from data provided by Tew 
(1986) and is labeled "SPDE-O." 
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Table 6-1 : Documented Engines 


Engine Type 

Operating Conditions: 

Reference 

Abbreviation 


Speed or frequency 
Mean pressure 
Working fluid 




Kinematic 

: Engines 


General Motors 
GPU3 

2500 rpm 
41 bar 
He 

Uriel! and Berchowitz 
(1984), pp. 37 and 39 

GPU3 

Mechanical 
Technologies, Inc. 
MODI (upgraded) 

4000 rpm 
150 bar 

h 2 

Richey (1986) 

MODI 

United Stirling AB 
P40 

4000 rpm 
150 bar 

h 2 

Tew (1983) 

pp. 80-86 

P40 

Technical 
University of 
Denmark , 
Stirling Total 
Energy System 

1500 rpm 
100 bar 
He 

Andersen (1979) 
pp. 50-52 

STES 

UK Consortium 
Stirling Engine 
( ot-conf i gur at ion) 

3000 rpm 
150 bar 
He 

Dunn et al. (1982) 

p. 68 

UK 

Ford Phillips 
4-215 

4000 rpm 
200 bar 

h 2 

Urieli and Berchowitz 
(1984), p. 30 

4-215 

General Motors 
4L23 

2000 rpm 
103 bar 

Martini (1982) 
p. 32 

4L23 



Table 6-1 : Documented Engines (continued) 



Operating Conditions: Reference 

Speed or frequency 
Mean pressure 
Working fluid 

Free Piston Engines 


Mechanical 

58 Hz 

Dochat (1985) 

Technologies, Inc. 

60 bar 


Engineering 

Model 

He 


Sunpower, Inc. 

60 Hz 

Berchowitz (1985) 

3 kW Generator 

25 bar 


Set 

air 


Sunpower, Inc. 

30 Hz 

Schreiber (1983) 

RE-1000 

70 bar 
He 


Mechanical 

trial parameters: 

Dochat (1985) 

Technologies, Inc. 

73 Hz 


Space Power 

75 bar 


Demonstrator 

He 


Engine 

extrapolated to design 

Dochat (1985) 

conditions: 
105 Hz 
150 bar 
He 


intended operation: 

105Hz 

150 bar 

He 

Tew (1986) 


Abbreviation 


EM 


Genset 


RE1000 


SPDE-T 


SPDE-D 


SPDE-0 
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7. RESULTS AND DISCUSSION 

The operating conditions of the heat exchangers in the Stirling engine 
data base are presented in terms of similarity parameters and conclusions 
are drawn about the conditions of fluid mechanics and heat transfer in 
heaters and coolers and in regenerators. The significance of various 
similarity parameters is discussed. Plots of heat exchanger operating 
conditions in terms of Re max , Re^ and Ar are used. 

7.1 Heaters and Coolers 

Entrance and exit losses . Table 7-1 gives a list of rough estimates of 
the points of the total heat exchanger pressure drop occurring at the heat 
exchanger entrance, in the duct of the heat exchanger, and at the exit. The 
estimates are based on the assumption of steady, unidirectional, turbulent 
flow in a rough pipe. Between 15$ and 50$ of the heat exchanger pressure 
drop occur at the pipe entrance and exit. Therefore entrance and exit 
effects are important. They could be estimated from an experiment such as 
the one by Taylor and Aghili (1984) and should be correlated in terms of the 
area ratios of contraction and expansion and Re max and Re^. 
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Engine 


Portion of the Pressure Drop 

Abbreviation 

(in percent) occurring... 

H - heater 

at the 

in the 

at the 

C - cooler 

entrance 

pipe 

exit 

GPU3 

H 

15 

68 

17 


C 

22 

53 

25 

MODI 

H 

14 

71 

15 


C 

14 

71 

15 

P40 

H 

15 

70 

15 


C 

16 

67 

17 

STES 

H 

22 

56 

22 


c 

24 

50 

26 

UK 

H 

7 

85 

8 


C 

24 

47 

29 

4-215 

H 

12 

75 

13 


C 

14 

71 

15 

4L23 

H 

14 

70 

16 1 


C 

13 

74 

13 

EM 

H 

16 

67 

17 


C 

16 

67 

17 

Genset 

H 

18 

62 

20 


C 

16 

62 

22 

RE1000 

H 

15 

68 

17 


C 

13 

72 

15 

SPDE 

H 

17 

65 

18 


C 

18 

62 

20 


Table 7-1 : Estimated pressure drop distribution 

in heaters and coolers. 

Estimated area ratios for the heater-regenerator and cooler-regenerator 
flow are listed in Table 7-2 to provide the similarity parameters for 
experiments on entrance and exit losses in oscillating flow. 
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Engine 

Abbreviation 

Heater-Regenerator 

Cooler-regenerator 

GPU3 

11 . 

ii . 

MODI 

21 . 

13 . 

P40 

40. 

17. 

STES 

32. 

16. 

UK 

29. 

5.9 

4-1 25 

30. 

18. 

4L23 

9.2 

18. 

EM 

13. 

9.4 

Genset 

7.6 

3.0 

RE1000 

7.8 

4.5 

SPDE 

13. 

9.0 


Table 7-2: Ratios of regenerator housing area to 

heater/cooler flow area. 

Developing profiles . Table 7~3 lists i/d values. For heaters and 
coolers with rectangular cross-section, the hydraulic diameter is used. In 
steady turbulent flow, the hydrodynamic entrance length is independent of 
Reynolds number and between 10 and 60 pipe diameters long (Incropera and 
DeWitt 1981, p. 377). Therefore, developing velocity profiles may be 
expected throughout many heaters and coolers. 
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Engine 

U/d) 

Abbreviation 

Heater 

Cooler 

GPU3 

81 .1 

42.7 

MODI 

96.9 

92.0 

P40 

93.7 

80.0 

STES 

51.0 

38.9 

UK 

219.0 

30.8 

11-215 

115.0 

96.7 

4L23 

88.6 

112.0 

EM 

81 .1) 

75.9 

Genset 

59.1 

50.9 

RE1000 

88.5 

77.6 

SPDE 

71 .0 

62.5 


Table 7-3 : Length to diameter ratios of heaters 

and coolers. 


Flow patterns and axial heat transfer . Figure 7-1 shows that the 
operating conditions of heaters and coolers fall within two decades of 
dimensionless frequency. The amplitude Reynolds numbers are all above the 
Reynolds number at which transition occurs in steady unidirectional flow. 
Most heaters and coolers operate in the turbulent flow regime but many are 
near transition. Therefore, it is important to understand the transition 
process and the characteristics of fluid mechanics and heat transfer in 
turbulent oscillating flow. Three heat exchangers operated under 
transitional conditions: the RE1000 cooler and the SPDE heater and cooler. 

The SPDE cooler is shifted from turbulent to transitional operation as the 
engine frequency is reduced. The same i3 expected if the Genset cooler or 
heater is operated at frequencies below the design point. Based on the 
discussion above, higher power dissipation and pressure drop than predicted 
by quasi-steady analysis are expected in the laminar and transitional 
regimes. Due to fluid inertia, there will be a phase shift between pressure 
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and mean velocity. Therefore, the instantaneous pressure drop cannot be 
determined from the instantaneous velocity and acceleration alone. 

The trial and design operating conditions of the SPDE fall on a 
straight line of slope 1 in the double logarithmic coordinates of Figures 
7-1 and 7-3. The ratio of the two similarity parameters is 

u max d 

Re max _ v _ it u max 

Re u>d 2 o)d 

0) t — 

Since d is constant and u max is proportional to u for any one engine, the 
ratio is constant. Therefore Re max is proportional to Re u if operating 
conditions are extrapolated assuming that only <d changes. 
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Figure 7-2 shows that the four heat exchangers in the transitional 
regime are operating with A R < 1 ; therefore, one might think that axial 
augmentation would be important. Estimates of this heat transfer 
augmentation for the SPDE heater and cooler indicate that the axial heat 
transfer is only about 0.1 - 1.0% of the heat transferred by the heat 
exchanger, however. This estimate was made using Watson's (1983) analysis 
and is based on the assumtion that the pipe walls are adiabatic. Further 
work is required to extend this analysis to non-adiabatic wall conditions as 
in the Stirling engine. 

Figures 7-2 and 7-M show that the relative amplitude of fluid 
displacement A R is the same for the trial and the design conditions, while 
the dimensionless frequency Re u is different. The trial and design 
conditions differ in the operating frequency and in the fluid density 
because of the different mean pressures. Therefore 



is affected while A R is not affected. 
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Compressibility effects. Table 7-4 shows that the Mach numbers are 


small In heaters and coolers, % ax < 0.075, with the exception of the UK 
heater (M max - 0.15). The density variation due to high velocity in the 
latter case is (Schlichting 1979, p. 10): 

f - \ Mmax 2 “ 0.012 

Therefore, high density variations due to high fluid velocity are not 
expected in the coolers and heaters of these Stirling engines. 


Engine 

Abbreviation 

Highest Mach 
Heater 

Number in the 
Cooler 

GPU3 

0.042 

0.034 

MODI 

0.067 

0.046 

P40 

0.064 

0.030 

STES 

0.035 

0.015 

UK 

0.153 

0.031 

4-215 

0.072 

0.043 

4L23 

0.027 

0.040 

EM 

0.020 

0.029 

Genset 

0.038 

0.026 

RE1000 

0.007 

0.020 

SPDE-T 

0.008 

0.008 

SPDE-D 

0.011 

0.009 


Table 7-4: Mach numbers in heaters and coolers. 
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7.2 Regenerators 


Figure 7~3 shows that the flow pattern in most regenerators may be 
expected to be complex. According to the observations by Dybbs et al. 
(1984), the flow would be turbulent or transitional for 8 of the 12 
regenerator operating points. Since their observations were made in steady 
state experiments, they may not hold in oscillating flow or may hold only 
for low Re^ and will probably only apply to matrices similar to those 
studied by Dybbs et al.(1984). 

Six of the 11 engines shown in Figure 7-4 have Ar < 1 which may limit 
the heat transfer from the matrix to the adjacent heater and cooler 
("holdup") . 

Mach numbers in the regenerators are M ^ y < 0.01 (Table 7-5). Choking 

of the regenerators is therefore not expected. At these low Mach numbers, 

the pressure drop in the regenerator matrices is expected to be independent 

of Mach number, according to the results of Beavers and Sparrow (1971) and 

% 

Benson and Baruah (1965) discussed in section 5.1. Since porosities of 
regenerators typically exceed 0.5 and Mach numbers are less than 0.2, 
choking in regenerator screens is not expected according to experimental 
results by Pinker and Herbert (1967, p. 16 and 17) which were discussed in 
Section 5.1 . 
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Engine 

Highest Mach Number 

Abbreviation 

in the Regenerator 

GPU3 

0.0078 

MODI 

0.0093 

P40 

0.0061 

STES 

0.0031 

UK 

0.0103 

4-215 

0.0055 

4L23 

0.0047 

EM 

0.0056 

Gen set 

0.0096 

RE 1000 

0.0038 

SPDE-T 

0.0017 

SPDE-D 

0.0024 


Table 7-5: Mach numbers in regenerators. 


7.3 Similarity Parameters 

After a review of the operating conditions of Stirling engine heat 
exchangers, several similarity parameters can be eliminated. Eckert numbers 
were calculated from 

Ec - (7 - D jj Mmax 2 

They range from 0.0066 (SPDE-D, heater) over 0.025 (MODI, heater) to 0.13 
(UK heater). Therefore, viscous dissipation is not considered to contribute 
significantly to the heat transfer, with the exception of the UK heater. 

The exact numerical value of the Eckert number is not important, provided 
that it is less than -0.1. 

A similar argument can be made for the Mach number , because Tables 7 4 

and 7-5 show that they are small, with the possible exception of the UK 

heater. Thus similarity with respect to Mach number is not required as long 

as M max < 0.1. Organ (1982a) challenged this argument with experimental 
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pressure drop results, but his data seem to support a strong frequency - 
(rpm-) dependence of pumping power. 

Pressure propagates through the heat exchangers fast compared to the 
cycle time. Table 7-6 provides estimates of the time required for 
propagation of pressure through the heat exchangers. It is assumed that 
pressure waves travel with the speed of sound which underestimates the 
propagation speed slightly. The data are given in degrees of crank angle. 
In all heat exchangers but the UK and 4-215 heater, pressure propagates 
through the Individual heat exchangers within 1/100 of the cycle time. 
Therefore, effects due to the finite propagation speed of pressure are not 
expected to be significant. 


Engine 

Heater 

GPU-3 

2.0 

MODI 

3.3 

P40 

2.9 

STES 

0.76 

UK 

5.2 

4-215 

4.5 

4L23 

2.0 

EM 

1.2 

Genset 

2.8 

RE1000 

1.2 

SPDE-T 

1.7 

SPDE-D 

2.4 


Cooler 

Regen. 

0.69 

0.24 

1.5 

0.83 

1.3 

0.51 

0.58 

0.094 

1.6 

0.64 

1.5 

0.43 

1.1 

0.16 

2.0 

0.86 

6.2 

3.1 

0.82 

0.52 

2.4 

0.54 

3.5 

0.78 


Table 7-6: Time estimates for pressure wave propagation 

through heat exchangers in degrees of crank angle. 
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8. CONCLUSIONS AND RECOMMENDATIONS 


In this report, important similarity parameters are introduced: 

Re max like Reynolds number in steady flow 

Re^, measure of unsteadiness 

a R» A Rh relative amplitude, important for heat transfer 

i/d, l h /d developing length effects 

Length ratios maintain geometric similarity. 

De like Re max for curved pipes 

Area ratios maintain similarity for entrance and exit losses. 

For porous media, the square-root of the permeability, /IT, may be 
the appropriate length scale to be used in Remax Re w 

Previous research results are discussed and 11 Stirling engine operating 
points are plotted in terms of the parameters. 

Due to the acceleration and deceleration of the fluid, velocity 
profiles and stability of laminar oscillating flow differ from steady 
laminar flow. Transition occurs at higher Re max as Re u is increased. 
Adequate turbulence models or correlations for turbulent oscillating flow do 
not exist. Results of past research on heat transfer are not generally 
applicable except the results for enhanced axial transport in oscillating 
laminar flow. It was shown that the heat exchangers of several engines are 
operating in or near the laminar-to-turbulent transition region while most 
heat exchangers operate in the turbulent regime. Several heat exchangers 
are long relative to the fluid displacement, A R < 1 . This means that some 
heat exchangers may add to the engine dead volume. 

The effect of compressibility of the working gas on pressure recovery 
can be accounted for by assuming an adiabatic expansion at the heat 


82 



exchanger and duct outlets. Choking is not expected due to the low Mach 
numbers encountered. Incipient shock formation was shown to be unlikely 
using a worst-case analysis. Pressure waves were calculated to pass through 
the heat exchangers quickly compared to the cycle time. 

A review of the literature has shown that immediate attention is needed 
in the following areas to assess the trade-offs between heat transfer and 

viscous dissipation in oscillating flow: 

1 . The process of transition must be understood more clearly. Work on 
this topic in tubular geometries will improve the understanding of 
heat transfer and pressure drop mechanisms in the heaters and 
coolers of Stirling engines. 

2. The means by which turbulence structure is affected by flow 
oscillation must be better understood, possibly leading to ways of 
modifying standard turbulence models so that turbulent oscillating 
flows can be accurately predicted. 

3. The effects of thermal and hydrodynamic entrance lengths on heat 
transfer and pressure drop in tubes and regenerator matrices must 
be better understood. 

H. A test of the permeability model should be made. It may provide a 
more general characterization of regenerator matrices than the 
characterization by hydraulic radius and porosity currently used. 

Once these topics have been addressed, non-sinusoidal oscillating flow 
should be studied experimentally to determine how general the results of 1-4 
are. A computational program can then use the experimental data from 1-3 to 
calculate heat transfer in oscillating flow with compression and expansion 
of the working gas. Its results should be confirmed experimentally. 
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APPENDIX A: A SURVEY OF PRESSURE, REYNOLDS-NUMBER AND MACH NUMBER 

VARIATIONS 

Figures of the pressure, Reynolds number and some Mach number 
variations are provided to give the reader the opportunity to assess the 
quality of the data base discussed in Chapter 7. Figures A1 and A2 show 
results of the Isothermal analysis of a GPU-3 (Youssef 1986). Figures A3 
through A1 4 were obtained by Schmidt analysis. 

- Flow is positive towards the cold end. 

- full line - hot end 

- dashed line » cold end 
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Figure A-l: Results of the isothermal analysis of a GPU 3: 
pressure and Reynolds nunbers 
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Figures A-3 through A-14: 


Schmidt Analysis Results 




Figure A-3: GPU 3 
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Figure A-4: MODI 
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Figure A-5 
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Figure A-6: STES 


99 





Figure A-7: UK 
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Figure A- 8: 4-215 
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Figure A- 9: 4L23 
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Figure A-ll: Genset 
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Figure A-12: RE1000 
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Figure A-14: SPDE-D 
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APPENDIX B: DERIVATION OF SIMILARITY PARAMETERS 

Normalization of the momentum equations 

To normalize the momentum equation, the following dimensionless 
variables are used: 


P - go 

p U* 


# p 


* U 

i m — 


- Ar 


Neglecting gravity, the momentum equation is: 




- 32 . , Ji r U 

P P 


+ Ml 5*. v *u* . _ + ^ u ^ 

t Q 3t^ L L p p 0 L 2 p 


klEo + Uk£fl - - £oHk ♦ v*V* 2 u* 

to^o mq Po P 


108 



For the oscillating 
length scale 
velocity scale 
time scale 


pipe flow problem, 
L - d/2 
U “ u m,max 

t Q - 1 /ta> 


the following scales are used: 
pipe radius 

amplitude of the mean velocity 
period of oscillation / 2ir 


^1! |ul + 1 H m.m ax d . _ 1 + V » v *2*» 

** v Q at 2 v 0 2 v Q p 


Re 

« at^ 2 


max u*. 


_# *►# 
V u 


- . v*7*»u* 


V*p* . 


where Re, 


max 


- d and Re - “ d ’ 


(0 




The dimensionless frequency Re w appears as a coefficient of the unsteady 
inertia term in the normalized Navier-Stokes equation and the amplitude 
Reynolds number Re max appears as a coefficient of the steady inertia term. 
Dividing by the maximum Reynolds number, the following form is obtained: 


-—a d |i £ + u** v*u* . 

2 u m,max 


# # 


-U- 


2 Vn 

Pm, max P 


v*7**u* 


Str + +* *+* . 

2 at^ 



In this form, the Strouhal number and the Reynolds number appear as the 
coefficients: 


Re max - U m > max d and Str 

max Vq 


to d 

u m,max 
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Normalization of the energy equation 


The energy equation Is i 


pc p ( i? + ^ ,VT) “ ft + ^' Vp + V ‘ (kVT) + * 


where * is the dissipation function, using Stokes’ hypothesis (e.g. 
Schlichting 1979, p.60): 


8u 


* " T U 3^7 " ^ 




The normalized variables are chosen as: 


t - u>t 


+# x 
x “ d/2 


-►# 

u - 


u 


u m f max 


p * . P ~ Pmp . a^ 
Po u m,max 


t* T - Tn 
‘ T h - T c 


T - d/2 V 


P “ P/Po 


k/k. 


P* - p/po 


'P' u po 


*■ - »• <!^> 2 


With these variables, the energy equation is: 

no 



Po c po P c p - T c ) -^r + (u mflDax u • ^2 ^ ~ T c )) 


3p* ■»# V* 

“ Po u m,max 2 "gpr + ^ u m r max u *’ ^72 p *^ p ° u m,njax 2 


* (Th - T c ) k Q ^ • (k* ^ T*) + po (Hnumax)* ** 


>' c p * [-^ |£* u-.vV] - g V max. j/2 |l£ + , W 2 J$*.yV 

^ u m , max c po^h * c ) c pO'^h T c ) 


Po c p d/2 u m>max 


V*-(k*V*T*) + HiMnax * 


Po c po d/2 (Tj^ ** T c ) 


*_ * \2 Re 3T* >#„#„# i 2 Re Ec 3p* . ># 

p C P LrT— “ at* u v T J " ~rP° af» + Ec u * 

ne max o 1 ne max dt 


* # 


V'p 


1 


Pr Re, 


V*-<h*V*T*) * S — ♦* 


max 


Re, 


max 


Finally, the normalized energy equation is: 


| | 

p*c p * (2 Re^ |£* + Re max u*-V*T*) - 2 Re^ Ec |£w ♦ Ec Re max S*-V*p* 


w 3t 


y*. (k*V*T*) » 

+ ♦ EC * 


Vi ^ 

where Ec - l Pj, max „ v is the Eckert number. 

c p (T h - T l ) 


Similarity parameters 1/d and Ar 

The similarity parameters are obtained from the normalization of the 
streamwlse component (x,u) of the momentum equation. Geometric similarity 
ensures similarity in the other components (y and z), e.g., the l/d - ratio 
in pipe flow ensures that the boundary layer growth in the hydrodynamic 
developing length is similar. 
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The set of parameters relevant for heat transfer in oscillating flow 
then consists of: */d, Re max , Re u , Pr, Ec and Ar. This set contains two 
parameters more than in steady flow. One of these parameters, however, is 
redundant for sinusoidal flow; if the velocity is 
u m - Um, max 3ln 

then the maximum displacement of an average fluid element is calculated by 
integration over half the cycle 


By 


2 x m>max - 0 J T/2u m,raax sin(a,t) dt 
the definition of Ar 


2 

u 


« „ 2 X m t maY _ 2 Um,max,_ 

A R A wt 

but Mmax . u m , max 1 1 . 4 u m . rna x_ 

DUt Re to d 2 to d 

0) 

therefore A R “ 2 I Re 

U) 

For other, non-sinusoidal fluid motion, the integration will yield a 
different proportionality factor than 1/2. In any case, Ar is a function of 
A/d, Re max and Re w . Therefore, the set of similarity parameters for heat 

transfer is reduced by one to 

Ar or ( fc/d) , Remax • R e u* ^* c an ^ * >r ‘ 

Experimental data will show whether Ar or A/d is more appropriate in 

correlating heat transfer results. 
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Calculation of the Reynolds number for regenerator matrices 

The hydraulic diameter for regenerator matrices In the Stirling engine 
data base Is calculated based on the definition used by Kays and London 
(198H, p. 8): 

dft/1 - $ A c /Ah 

where A c - flow cross-sectional area 
Ah - total heat transfer area 
£ - heat exchanger length 

d^ - x void volume/heat transfer area 
$ - porosity 

Figure B-1 shows a cross-section of two layers of a woven-screen 
regenerator, where s is the wire spacing (i.e., the Inverse of the screen 
pitch) and d w is the wire diameter. 


s 





CROSS-SECTION FRONT VIEW (INFLOW DIRECTION) 

Figure B-1 : Schematic of a wire mesh regenerator screen. 


To a first approximation, the control volume (c.v.) contains a solid 
volume 
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The porosity Is defined as the ratio of void volume V v to total volume V t . 



The heat transfer area in the control volume is 
Ah - * d w s 

Then the hydraulic diameter is (Kays and London 1984, p. 8) 


d h - 


* Lt - * ±Jk . 

♦ ♦ *h 


♦ la 

♦ " 1 Ah 


d w 




This is the expression commonly used for the hydraulic diameter. 

The characteristic velocity in the Reynolds number is calculated as the 
average velocity in a cross-section with effective flow area 

A c “ ♦ s* 

Then the characteristic velocity is 


V 

u m - * 

where V is the superficial mean velocity, i.e., the mean velocity of the 
flow just before entering or after leaving the matrix. 
u mA c - VAf 

where Af is the total regenerator frontal area, i.e., the cross-sectional 
flow area when the matrix is removed. 
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APPENDIX C: VELOCITY PROFILES IN LAMINAR FLOW 


The plots below show velocity profiles for laminar, incompressible, 
fully developed, oscillating flow for eight dimensionless frequencies. The 
velocity profiles were calculated from the equations given by Uchida (1956) 
for pulsatile flow by taking the limit as the mean velocity goes to zero. 
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Figure C-l: Velocity 


NORIIflUZED VELOCITY NORtWLIZED VELOCITY 
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APPENDIX D: EQUATIONS OF OBSERVATIONS OF TRANSITION 


This appendix provides the equations describing the lines in Figure 3-8. 
The ranges of applicability are based on the range over which the original 
references provide experimental data to support these equations. 


Reference 

Equation 

Range 

Grassmann and Tuma (1979) 

Remax - W Re*, 0 * 75 

42 i Re u £ 520 


R e max “ 15300 

520 < Re u S 943 

Iguchi et al. (1982) 



laminar ■* transitional 

R^max * ^00/Re w 

100 < Re w < 800 

transitional -*• turbulent 

^ e max “ 800/Re w 

100 < Re u < 800 

Ohmi et al. (1982) 

Re max “ 800/Re w 

18.4 S Re u S 546 

i 

Park and Baird (1970) 

i 

R e max “ 1 

35 S Re u $ 1000 

Sergeev (1966) 

R e max * 700/Re^ 

16 S Re u S 1600 


118 



APPENDIX E: EFFECTS OF PRESSURE PROPAGATION 

Resonance 

A necessary condition for resonance in Stirling engines is that the 
flow length between pistons is greater than half the wavelength (X) of 
pressure propagation. 



where a - speed of sound 

f - engine frequency 

The smallest X is obtained for small a (low temperature) and large f. For 
the only air-charged engine in the data base (Genset), the smallest X is: 


7rT~ _ /I. HI x 286.7 X 367 m/sec _ 6> „ 2m 
f “ 60 sec" 1 


For He and H 2 -eharged engines the worst case would be the He-charged SPDE-D, 
the smallest X is: 

/YRT ~ /I. 66 x 2079 x 350 m/sec _ 10<5m 
* " f " 105 sec” 1 

The flow length between pistons may be found in the literature for the HL23 
(Martini 1982, p. 32) as 1.4m. It can be estimated for the GPU3 from Uriel! 
and Berchowitz (198H, p. 36-37) as - O.Hm and for the UK engine from Dunn et 
al. (1982, p. 68) as > 0.78m. For all other engines in the data base it is 
expected to be less than 3m, estimated from the exterior dimensions of the 


engines. 
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The half-wavelengths, A/2, for both cases exceed the flow length 
between pistons, therefore, resonance is not expected. 

Shock incipience 

A shock forms when pressure waves interfere constructively. This may 
be predicted by the method of characteristics. Consider a piston that 
oscillates sinusoidally in a tube (Fig. E-1). 



Figure E-1 : Oscillating piston in a tube 

As the piston moves to the right, it continuously sends pressure waves to 
the open end. These can be represented in the t-x-plane as shown in Fig. 
E-2. Two pressure waves that were emitted by the piston at different times 
may interfere if the one emitted later travels faster (smaller slope of the 
t vs. x line). 
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The lines describing the pressure wave propagation are characteristics. If 
the piston velocity is small compared to the velocity of sound, the slopes 
of both lines will be nearly equal: 


(— ) - 

v dx ii Up + a 

,dt. ,dt. 

If a > Up, then (^) j - (^ n 

Then the two lines will intersect far from the piston. 

In an engine, the maximum piston velocity is smaller than the fluid 
velocity in the heat exchangers because the flow area is less than the 
piston area. Therefore, a worst case is obtained if the piston velocity is 
set to the highest fluid velocity found in the data base. The 
characteristics in Fig. E-2 are most likely to intersect if the period of 
the sinusoidal piston motion is small, i.e.. the worst case should have the 
highest frequency. Due to the temperature gradient in an engine, the speed 
of sound varies. For the present simplistic analysis, the temperature is 
assumed to vary linearly with x* 

The analysis is one-dimensional and ignores contractions, expansions, 
bends and regenerator matrices in the flow path. A fictitious worst case is 
composed of all the most favorable conditions for shock incipience 
encountered in the data base: 


piston 

velocity amplitude 

Mmax ” 0.15 

UK-heater 

frequency of oscillation 

to • 1 05 Hj 

SPDE-D 

engine 

flow length 

L - 1.5m 




L > 1.39m 

*«L23 
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The temperature was assumed to vary linearly with x from 300K to 1050K. 
Consequently, the speed of sound varies as: 


a(x) - /YRT(x) 

- ZyR/L [(T h " T c )x + T c ] 

Then the two characteristics are described by 

t# I “ 4F Th - T c t /(T h “ T 0 )x* + T c L - ^L] 
t # n " t 0 - Th \ Tc [/(Th ** T c )x ♦ T C L 

♦ U PQ^" l 1 ' ln(/(T h - T c )x * T C L ♦ u po )}]** 

where u p0 is the piston velocity at t 0 , when the pressure wave leaves the 
piston and * denotes shock incipience due to interference of pressure waves. 

Test cases for air, hydrogen and helium were calculated for different 
values of t 0 but shock incipience was not observed in any of these test 
cases (e.g., Fig. E-3). 
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APPENDIX F: AUGMENTATION OF AXIAL TRANSPORT 


From Watson (1983, p. 241), the augmentation coeffficient is 
Au - ¥ (Wo, Pr) Vf* Wo - /Re u 


with 

¥ 


1 B, (Wo) 

" B a (Wo /Pr) 

8ti z ( 1 - Pr* /z ) B,(Wo) 


where 


B a (Wo) 


Wo B w (Wo) ♦ B* (Wo) 
B' (Wo) 


and 


_ x Wo*B(Wo) + B* (Wo) - WoB"(Wo) - Wo a B"»(Wo) 
B * lWo; " Wo^B' (Wo) 

B n - (ber*(Wo) + bei*(Wo)) 


Ber and bei are the Kelvin functions. A series solution for B and its 
derivations is provided by Watson (1983. P« 241). 

Based on this series solution, the plot of the augmentation function ¥ 
shown below was generated (Figure FI). The axial heat transfer can be found 
as discussed in section 4.2 if an axial temperature gradient is assumed. 
Sample calculations are not included here because they involve proprietary 
geometry data. 
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